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A  MODIFIED  PRONY  METHOD  APPROACH  TO  ECHO-REDUCTION 
MEASUREMENTS  OF  TIME-LIMITED  TRANSIENT  SIGNALS 

INTRODUCTION 

The  current  procedure  for  making  echo-reduction  measurements  of 
acoustic  panels  at  the  Underwater  Sound  Reference  Detachment  of  the  Naval 
Research  Laboratory  employs  steady-state  signal  processing  of  measurements 
made  in  the  anechoic  tank.  The  anechoic  tank  allows  measurements  to  be 
performed  at  various  temperatures  and  hydrostatic  pressures.  However,  it 
restricts  panel  width  to  76.2  cm  due  to  the  size  of  the  access  port. 


Fig.  1  -  Panel  measurement  geometry 


Figure  1  illustrates  the  geometry  of  the  measurements.  The  panel  is 
insonified  by  a  projector  170  cm  in  front  of  the  panel,  and  the  incident 
and  reflected  signals  are  monitored  by  a  probe  34  cm  in  front  of  the  panel. 
For  a  step  sinusoid  signal  from  the  projector,  the  output  of  the  probe, 
illustrated  in  Fig.  2,  consists  of  426-us  segment  of  incident  signal 
followed  by  a  141-us  segment  of  both  incident  and  reflected  signals 
followed  by  the  arrival  of  the  diffracted  signal  from  the  panel  edges. 

In  order  to  eliminate  the  undesired  diffracted  signal,  the  measurements 
are  time  limited  to  approximately  140-ys  for  the  reflected  signal.  This 
time  limit  and  the  constraint  that  the  signal  reach  its  steady-state  re¬ 
sponse  currently  restrict  measurements  to  frequencies  above  15  kHz.  To 
obtain  measurements  below  15  kHz  it  must  be  possible  to  extrapolate  the 
steady-state  response  from  the  transient  portion  of  the  signal. 
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Fig-  2  -  Time  segment  output  of  probe  for  step  sinusoid  signal 


The  Prony  [1]  method  first  published  in  1795  allows  such  an  extra¬ 
polation  and  has  been  investigated  as  a  means  of  allowing  low-frequency 
panel  measurements.  In  this  report  a  modified  Prony  method  is  presented, 
which  allows  measurements  to  be  performed  down  to  2  kHz,  along  with  the 
results  of  measurements  on  various  panels  to  illustrate  the  ability  of 
the  modified  approach. 

THEORY 

Any  signal  that  can  be  represented  by  a  function  that  is  the  solution 
of  a  set  of  constant-coefficient  linear  differential  equations  has  an  ex¬ 
pansion  given  by  [2] 

N 

6U)  »  £  A  -  expU  ~t)  (1) 

i-i  3  3 

where  W  is  the  order  of  expansion  and  Ay  and  4  j  sure  respectively  the 
amplitudes  and  poles  to  be  determined.  Since  fait)  is  real,  both  Ay  and 
-6 j  must  be  real  or  appear  in  complex  conjugate  pairs.  The  real  part  of 
the  pole  -6y  is  the  time  constant  associated  with  the  transient  response, 
and  the  imaginary  part  is  the  angular  frequency  of  the  j th  component. 

The  phase  information  is  contained  in  the  amplitudes  Ay. 

There  are  two  immediate  difficulties  in  obtaining  the  expansion  in 
Eq.  (1) .  First  there  is  the  problem  of  the  nonlinearity  of  the  expan¬ 
sion  in  -&j  and  second,  there  is  the  indeterminacy  of  the  order  of  ex¬ 
pansion  W.  The  first  problem  was  solved  by  Prony  and  his  method  will 
be  described  below.  The  indeterminacy  of  the  order  of  expansion  will 
be  discussed  after  the  basic  equations  have  been  derived. 

Prony  found  that  the  desired  expansion  could  be  obtained  by  consid¬ 
ering  a  uniformly  sampled  version  of  Equation  (1)  then  becomes 

N 

6(kA)  =  7  A.  exp(4;fe A)  ,  (2) 

y=i  J  J 

where  d  is  the  data  sampling  time  interval  and  k  is  an  indexing  integer 
running  from  0  to  M  such  that  Aid  is  equal  to  the  observation  time. 


Let 


(3) 


Zj  =  expU-L)- 
N  fe 

Then  (feA)  =  £  A -z  .  .  (4) 

y=i  3  3 


Equation  (4)  does  not  represent  a  linear  set  of  equations  in  zy.  However, 
Prony  noted  that  solutions  to  zy  could  be  found  by  first  defining  an  W-C/i 
order  polynomial  in  z  whose  roots  are  zy.  That  is 

N  -  U 

P(z)  =  l  CL-Z  =  .n  (z-z  .)  ,  (5) 

i=o  -*’1  3 


where  =  1.  The  ctj_  coefficients  have  one  more  degree  of  freedom  than 
the  Zj;  hence  the  arbitrary  setting  of  equal  to  unity.  Using  Eqs.  (4) 
and  (5)  to  evaluate  the  following  expression: 


W 

l  CL L  4lU  +  £)A ], 
l~  0 


For  £  =  0  to  M-N  yields 


W  W  /  {(  ;+/\ 

X  ^UU  +  £)*]  =  .laA  X  Ayzy  ) 

-t=0  -1  x.=0  \i=i  3  3  / 

■aa4JIv/) 


v/  P|2y'  ■  «• 


Or 


N-1 

l  a  .  rf[(£  +  -£■)  A]  =  C(N  +  £)A] 
£=0 


(6) 


(7) 


since  =  1.  Equation  (7)  yields  M-W+ 1  linear  equations  in  the  W  unknown 
polynomial  coefficients.  For  M=ZN-1  the  equations  can  be  solved  to  obtain 
the  oi£  which  completely  define  the  polynomial  in  Eq.  (5) .  The  roots  of  the 
polynomial  can  be  found  to  provide  the  required  zy.  These  zy  are  then  sub¬ 
stituted  into  Eq.  (4)  to  determine  the  amplitudes  Ay.  The  expansion  is 
then  completed  by  transforming  the  zy  back  into  the  complex  4  plane  using 


4y  *  A-1  In  Zy.  (8) 

Prony' s  method  isolates  the  nonlinearity  of  the  expansion  to  Eq.  (5). 
However,  the  method  introduces  a  third  difficulty  through  Eq.  (8). 


Negative  real  Zj  cannot  be  transformed  back  into  the  complex  4  plane. 
Although  these  represent  nonphysical  terms,  noise  contaminated  data  will 
occasionally  produce  negative  real  zj.  When  this  occurs,  one  must  be 
satisfied  with  the  expansion  in  the  form  of  Eq.  (4) . 

In  addition  to  negative  real  Zj,  nonphysical  roots  occur  as  a  conse¬ 
quence  of  the  indeterminacy  of  the  order  of  expansion.  If  fewer  poles  are 
specified  in  the  expansion  than  actually  exist  in  the  data,  then  the  poles 
returned  by  the  algorithm  deviate  from  the  true  poles  and,  in  most  cases, 
the  true  poles  will  not  be  returned  at  all.  If  more  poles  are  specified 
than  actually  exist,  the  algorithm  will  return  extraneous  poles  in  addi¬ 
tion  to  a  set  that  deviates  from  the  true  poles.  The  presence  of  the 
extraneous  poles  adversely  affects  the  calculated  amplitudes  of  the  true 
poles.  Attempts  to  systematically  determine  the  order  of  expansion  by 
checking  for  linear  dependent  columns  in  the  matrix  generated  by  the 
Prony  difference  equations  or  by  analyzing  the  eigenvalues  of  the  matrix 
have  failed  in  the  presence  of  noise  [3] . 


In  the  modified  method  to  be  developed  in  the  remainder  of  this  section, 
the  problems  of  negative  real  zj  and  the  indeterminacy  of  the  order  of  ex¬ 
pansion  will  be  circumvented.  This  will  be  accomplished  through  the  intro¬ 
duction  of  a  priori  poles  and  will  be  discussed  after  the  equations  are 
developed.  However,  first  a  modified  least  square  extension  of  the  prony 
method  will  be  developed. 


The  solution  of  Eq.  (7)  requires  at  least  2 N  sampled  data  points.  For 
the  case  where  M  =  ZN-1,  the  solution  will  pass  directly  through  the  sampled 
data  points  since  the  problem  is  exactly  specified.  However,  in  general, 
the  number  of  available  data  points  exceeds  ZN  and  a  type  of  least  square 
solution  is  used.  This  is  most  conveniently  done  by  performing  a  pseudo¬ 
inverse  of  Eqs.  (7)  and  (4).  Equation  (7)  in  matrix  notation  is 


Ax  •  y, 


where 


(9) 


^0  ^1 
h  h 


h 

h 
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*N-l 

*Wl 


l  ^M-N  ^M-N+l  ^M-N+2 


*M-1 
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Following  the  same  procedure  the  least  square  extension  of  Ea.  (4)  is 


N 


z/ilkAi 


(12) 


for  i  5  1  to  M. 


This  least  square  extension  of  the  Prony  method  was  examined  by 
McDonough  [4] ,  who  observed  that  the  normal  equations  become  increasingly 
ill  conditioned  as  the  sampling  rate  increases.  The  source  of  the  diffi¬ 
culty  lies  in  the  data  interval  spanned  by  each  difference  equation.  The 
interval  spanned  is  NA;  and  as  the  sampling  rate  1/A  increases,  the  spanned 
interval  decreases.  As  smaller  portions  of  the  signal  are  used  in  each 
difference  equation,  the  noise  in  the  signal  has  a  greater  deleterious 
effect  on  the  pole  estimates.  In  order  to  eliminate  this  effect  Beatty 
[5]  uses  every  pth  sampled  data  point  instead  of  adjacent  sauries  to 
satisfy  the  difference  equations.  This  results  in  a  data  span  of  Np A 
for  each  difference  equation.  Then,  as  the  sampling  rate  is  increased, 
the  value  of  p  is  increased  to  maximize  the  interval  spanned.  All  of  the 
data  are  used  since  each  difference  equation's  initial  data  sample  is 
l$(feA)  where  k  is  an  incrementally  stepped  integer.  For  this  modification 
the  equations  are  developed  below. 


Let 


N 

( feA  +  npA)  =  l  A  -  exp[A-(fe  A  +  up  A  J  ] 
j*  l  J  J 


(13) 


where  k,  n,  and  p  are  all  positive  integers  and  p  is  the  introduced  data 
increment.  Let 


Z  y  =  eXp(AyPA)  (14) 

and 

{ 5!fc,n)  =  rfffeA  +  npA j .  (15) 

Then  Eq.  (13)  becomes 

N  W 

=  J  Ay  eXpUyfeA)  Zyn  =  J  B.z'.n  ,  (16) 

j  ~  1  j  —  1/J 

where  By  s  Ay  zxplAjk A),  since  the  difference  equations  do  not  depend 
upon  By  (fe  is  a  constant  in  each  difference  equation) ,  the  difference 
equations  are 

N-l 

l  Hk,n)a  *  -6(M>  (17) 

n= o  1 
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for  k  =  0  to  M-W.  The  least  square  extension  of  Eq.  (17)  is  obtained  by 
forming  the  pseudo-inverse  resulting  in 

N-l  1M-N  )  M-N 

I  «  \  l  i Uk,n)Uk,l)\*  -  l  Uk,N)&ik,l)  (18) 

n=o  n  { k=o  ’  k=o 

for  i.  -  0  to  SI-1.  The  solution  of  Eq.  (18)  yields  the  coefficients  an 
which  define  P(z').  The  roots  of  P  ( z"  )  are  Z 'j.  The  Z may  be  transformed 
back  into  the  complex  4  plane  using 


4  •  =  (pA)'1  In  Zy.  (19) 

Once  the  4y  cure  computed  they  may  be  substituted  into  Eq.  (3),  and  then 
Eq.  (12)  can  be  used  to  obtain  the  amplitudes  Ay.  As  before,  the  pro¬ 
cedure  fails  if  a  negative  real  Z'j  is  obtained;  except  in  this  case  not 
even  the  amplitudes  can  be  obtained. 

In  the  case  of  echo-reduction  measurements  the  useful  information  is 
contained  in  the  amplitude  of  the  steady-state  driving  frequency  component . 
Under  this  circumstance  the  Prony  method  can  be  modified  by  the  introduc¬ 
tion  of  a  priori  poles  to  circumvent  the  problems  of  negative  real  z'j  and 
the  indeterminacy  of  the  order  of  expansion.  If  the  algorithm  is  con¬ 
strained  to  find  the  correct  a  priori  poles,  the  remaining  poles  used  in 
the  expansion  are  used  as  curve-fitting  poles.  There  is  no  requirement 
that  the  curve-fitting  poles  have  any  physical  significance.  Then  the 
only  requirement  on  the  order  of  expansion  is  that  there  be  a  sufficient 
number  of  curve-fitting  poles  such  that  the  mean  square  deviation  between 
the  waveform  and  it's  expansion  be  below  some  arbitrary  small  value. 

The  problem  of  negative  real  z'j  is  eliminated  by  using  an  odd  integer 
for  the  data  increment  p.  Then  Eq.  (12)  may  be  used  directly  to  obtain 
the  amplitudes  since 


and  negative  real  Zy  are  handled  with 


Z  • 

J 


The  modified  equations  are  developed  below. 


Let 


(21) 


HU'  ) 


(22) 
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be  the  polynomial  generated  by  the  R  a  priori  poles  where  =  1 .  since 
N(z')  must  be  a  factor  of  Prony's  polynomial,  P(z')  in  Eq.  (5)  can  be 
written  as 

N-R 

P(z')  =  N(z')  l  Bz'n  , 
n=o  n 

where  =  1-  Writing  Eq.  (23)  explicitly  yields 

N-R 


l  V'*  -  z  V  l  tnz'n . 

l=o  *■  1=0  n=o 


Equating  similar  powers  of  z'  yields 


where  =  Bn-1  =  •  •  •  =  0A/-R+1  =  0  an*3  3-1  =  3-2  = 
Substitution  of  Eq.  (25)  into  Eq.  (-17)  yields 


(25) 


=  B-R  =  0. 


iV-1 

1 

4.=0 


l  *(M)  j  I  bzh-i]m 


(26 ) 


which  may  be  rewitten  as 

W-R-l 

■i 

Let 


■  R- 1  (  R  ]  r 

i  8r<  l  bMk,l+l)\  =  -  l  bMk,N-K+j) 
1=0  ^\z=o  4  )  iso  j 


(27) 


F(fe,n)  «  [  6,rflfe,n+/) 


(28) 


Then  Eq.  (27)  becomes 

N-R-l 

J  P(fe,X)  =  -F  { fe,  W-  R ) .  (29) 

Z=0 


Forming  the  pseudo-inverse  results  in 

W-R-l  (M-N  ) 

X  6  -s  l  F{k,Z)F(k,j}\  = 

i*o  Hk=o  ) 


M-N 

I  F(fe,N-R)F(fe,yj 
k=0 


(30) 
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for  /  =  1  to  W-R-l.  Equation  (30)  is  solved  to  obtain  the  coefficients  B^. 
These  6^  are  then  used  to  generate  the  reduced  Prony  polynomial  whose  roots 
z'-  are  the  curve-fitting  poles.  The  curve-fitting  poles  together  with  the 
aJpriori  poles  are  transformed  using  Eqs.  (20)  and  (21).  The  computed  Zy 
are  then  used  in  Eq.  ( 12)  to  find  the  amplitudes  Ay.  The  amplitude  of  tne 
steady-state  driving  frequency  pole  is  then  used  to  compute  the  echo 
reduction  in  a  manner  to  be  discussed  later. 

COMPARISON  OF  MODIFICATIONS 

In  order  to  test  the  effectiveness  of  the  various  modifications  on  the 
type  of  signal  to  be  encountered  in  panel  measurements,  the  waveform  in  Fig. 

3  was  generated.  The  waveform  simulates  the  reflection  of  a  3-kHz  step 
sinusoid  from  a  0.95-cm  thick  infinite  steel  plate.  It  was  computer  gener¬ 
ated  by  successively  adding,  with  suitable  time  delays,  the  multiple  internal 
reflections  that  are  transmitted  back  through  the  face  of  the  plate.  Seven 
data  files  were  constructed  by  sampling  the  waveform  at  1  MHz  and  adding 
various  levels  of  random  noise.  The  first  200  ps  of  each  data  file  was 
then  analyzed  by  the  Prony  method  in  six  different  manners. 
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Fig.  3  -  Simulated  reflection  of  a  3-kHz  step  sinusoid 
from  a  0.95-cm  thick  infinite  steel  plate 

The  six  expansions  performed  were  divided  into  two  sets  of  three,  one 
of  which  used  three  poles  in  the  expansion  while  the  other  used  fifteen  poles 
Three  pole  expansions  were  used  since  the  waveform  has  a  known  three-pole 
expansion  consisting  of  a  complex  conjugate  pair  representing  the  steady- 
state  driving  frequency  and  a  real  pole  associated  with  the  transient 


response  of  the  plate-  The  choice  of  fifteen  poles  was  arbitrary.  Each 
set  contained  one  expansion  with  no  a  priori  poles  and  a  data  increment  of 
one  one  expansion  with  no  a  priori  poles  and  a  data  increment  of  eleven, 
and  an  expansion  with  two  a  priori  poles  and  a  data  increment  of  eleven. 

The  a  priori  poles  entered  were  the  steady-state  driving  frequency  poles 
and  since  200  data  points  were  used,  all  expansions  used  least  square 
methods. 

The  amplitude  of  the  steady-state  driving  frequency  poles  or  the  poles 
closest  to  the  driving  frequency  when  no  a  priori  information  was  entered 
were  used  as  a  measure  of  the  accuracy  of  the  expansion.  The  results  are 
plotted  in  Fig.  4  where  the  correct  amplitude  is  0.426  and  the  two  expan¬ 
sions  with  a  data  increment  of  one  were  not  plotted  since  the  expansions 
failed  in  most  cases  to  obtain  any  poles  close  to  the  driving  frequency. 
The  results  indicate  that  unless  one  has  a  high  signal  to  noise  ratio  the 
only  method  that  obtains  useful  information  is  the  use  of  both  a  priori  and 
curve  fitting  poles.  In  general  the  more  a  priori  information  supplied  and 
the  more  curve  fitting  poles  used  the  better  the  results. 
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Fig.  4  -  Amplitudes  of  3-kHz  component  of  simulated  panel  reflection 
obtained  by  Prony  algorithm  with  various  modifications 
and  noise  levels.  Correct  amplitude  is  0.426. 


SELECTION  OF  INPUT  VARIABLES 


The  modified  Prony  algorithm  contains  three  user  supplied  variables: 
the  data  increment,  order  of  expansion,  and  length  of  time  window. 
Unfortunately,  no  strict  limits  on  the  variables  can  be  set  since  they 
depend  on  the  complexity  of  the  waveform  being  analyzed  and  the  signal- 
to-noise  ratio.  However,  given  the  constraints  under  which  the  method 
will  be  used  in  making  echo-reduction  measurements,  several  useful  comments 
can  be  made.  In  this  section  the  required  values  of  the  variables  for 
echo-reduction  measurements  will  be  investigated. 

Data  Increment 

In  order  to  investigate  the  variables  with  the  actual  signals  to  be 
encountered  in  panel  measurements,  a  set  of  waveforms  was  obtained  for 
the  reflection  of  a  step  sinusoid  from  a  1. 27-cm- thick,  76-cm-square 
steel  panel.  The  waveforms  were  obtained  by  placing  a  probe  5  cm  and  a 
projector  170  cm  in  front  of  the  steel  panel.  A  1-ms  pulse  from  the 
projector  produced  a  6 7- us  segment  of  direct  signal,  at  the  probe, 
followed  by  250  ps  of  incident  plus  reflected  signals  before  the  arrival 
of  the  diffracted  signal  from  the  panel  edges.  The  output  voltage  of 
the  probe  was  sampled  at  1  MHz,  and  100  separate  recordings  were  averaged 
to  reduce  the  incoherent  noise  level.  A  second  set  of  measurements  were 
made  without  the  steel  panel  in  place  to  obtain  a  long  recording  of  the 
incident  signal.  The  two  waveforms  were  then  directly  subtracted  to 
yield  the  reflected  waveform  from  the  steel  plate.  Figure  (5)  illustrates 
the  waveform  of  the  steel  panel  reflection  at  3  kHz. 
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Waveform  obtained  by  the  reflection  of  a  3-kHz 
step  sinusoid  from  a  1.27-cm-thick  steel  plate 


The  modified  Prony  method  was  used  to  analyze  the  3-kHz  waveform  to 
obtain  the  amplitude  of  the  steady-state  driving  frequency  poles  that 
were  entered  a  priori.  A  15-pole  expansion  was  used,  and  the  time  window 
was  reduced  in  10-ys  steps  from  250  to  110  ys.  Four  different  data  in¬ 
crements  were  used,  and  the  results  are  illustrated  in  Fig.  (6) .  Values 
for  the  three  largest  data  increments  do  not  span  the  entire  time  scale. 
This  is  due  to  the  requirement  of  a  minimum  number  of  data  points  for  the 
expansion  as  determined  by  the  order  of  expansion  and  the  data  increment. 
The  minimum  number  of  data  points  required  by  the  algorithm  is  given  by 

Min.  #  pts.  =  (Data  Increment  +1)  *0rder  of  Expansion.  (31 


With  a  15-pole  expansion  the  minimum  number  of  data  points  for  a  data 
increment  of  11  is  180.  Since  the  waveform  was  sampled  at  1  MHz,  the 
minimum  time  window  for  a  data  increment  of  11  is  180  ys. 


TIME  WINDOW  (ys) 


Fig.  6  -  Amplitude  of  the  steady-state  driving 

frequency  poles  for  waveform  in  Fig.  (5) 
obtained  with  15-pole  expansions  and 
various  data  increments 

The  results  in  Fig.  (6)  are  consistent  for  large  time  windows  but  vary 
as  the  time  window  is  reduced.  Since  in  general  the  largest  possible  data 
increment  should  be  used,  the  appropriate  portions  of  Fig.  (6)  have  been 
reproduced  in  Fig.  (7). 
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Fig.  (7)  -  Segments  of  data  in  Fig.  (6)  illustrating 
effect  of  minimum  time  window 

Figure  (7)  illustrates  the  one  exception  to  the  largest  possible  data 
increment  rule.  The  values  at  180  us  for  a  data  increment  of  11  and  at 
150  ys  for  a  data  increment  of  9  show  marked  deviations  from  the  average 
value.  Both  of  these  cases  correspond  to  minimum  values  for  the  time 
window  of  the  associated  data  increment.  The  difficulty  lies  in  the  matrix 
generated  from  the  Prony  difference  equations.  If  no  a  priori  poles  are 
used,  the  minimum  number  of  points  in  Eq.  (31)  generates  an  n*n  matrix 
where  n  is  the  order  of  expansion.  Since  this  is  a  square  matrix,  no  least 
square  technique  is  required  to  solve  the  equations.  However,  if  a  priori 
poles  are  used,  the  matrix  generated  is  an  (n-.'ljxk  matrix  where  K  is  the 
number  of  a  priori  poles  used  in  the  expansion.  The  matrix  represents  an 
overdetermined  set  of  equations  for  n-A.  unknown  coefficients,  and  a  least- 
square  technique  must  be  used  to  solve  the  equations.  When  K  is  small  in 
comparison  to  n  or  when  the  matrix  is  nearly  a  square  matrix,  the  least- 
square  technique  introduces  considerable  error  into  the  calculation.  All 
of  the  measurements  in  Figs.  (6)  and  (7)  were  done  with  an  order  of 
expansion  of  15  and  with  two  a  priori  poles.  By  reducing  the  value  of 
the  data  increment  when  the  minimum  time  window  is  approached,  the  matrix 
is  no  longer  nearly  a  square  matrix,  and  the  least-square  technique  returns 
consistent  values  as  illustrated  in  Fig.  (7).  Thus  the  largest  data 
increment  should  be  used  except  when  the  minimum  time  window  is  approached, 
and  then  the  next  lower  value  should  be  used. 


Figure  (7)  also  indicates  that  for  time  windows  of  less  than  150  us 
an  insufficient  portion  of  the  waveform  is  being  used  to  yield  useful 
results. 
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Order  of  Expansion 


The  order  of  expansion  is  the  most  difficult  variable  to  determine. 
There  must  be  a  least  as  many  poles  in  the  expansion  as  there  are  in 
the  signal.  However,  a  knowledge  of  the  structure  of  the  signal  does 
not  guarantee  correct  results.  As  illustrated  in  Fig.  (4)  a  3-pole 
expansion  with  two  poles  entered  a  priori  was  sufficient  for  a  waveform 
that  was  known  to  have  only  three  poles  when  the  signal-to-noise  ratio 
was  50  dB.  When  the  signal-to-noise  ratio  was  lowered,  the  3-pole  ex¬ 
pansion  yielded  incorrect  results. 

In  general  the  more  poles  used  in  the  expansion  together  with  the 
a  priori  poles,  the  better  the  results.  Unfortunately  the  larger  the 
order  of  expansion,  the  longer  the  running  time  for  the  program.  The 
optimum  value  is  strongly  dependent  on  the  complexity  of  the  waveform, 
the  signal-to-noise  ratio,  and  the  length  of  the  time  window.  To  get 
some  idea  of  the  required  order  of  expansion  for  echo-reduction 
measurements,  the  waveform  in  Fig.  (5)  was  again  analyzed.  Five  differ¬ 
ent  orders  of  expansion  were  used  on  time  windows  that  varied  from  250 
to  150  ys.  The  data  increment  was  determined  by  the  results  of  the 
previous  section. 


In  Fig.  (8)  the  amplitude  of  the  steady-state  driving  frequency  pole, 
which  was  entered  a  priori,  has  been  plotted  against  the  length  of  the 
time  window.  The  15-pole  expansion  deviates  by  less  than  0.1  dB  over 
the  entire  time  window  span.  The  12-pole  expansion  is  consistent  down 
to  a  time  window  of  180  ys.  However,  the  remaining  expansions  yield 
inconsistent  results  and  vary  from  one  time  window  to  the  next.  This 
indicates  that  the  order  of  expansion  must  be  at  least  15  for  echo- 
reduction  measurements  of  simple  homogeneous  plates  and  may  have  to  be 
higher  for  nonhomogeneous  plates. 
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Amplitude  of  the  steady-state  driving  frequency  poles  for  wave¬ 
form  in  Fig.  (5)  obtained  with  various  orders  of  expansion 
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Time  Window 


The  minimum  time  window  necessary  for  the  algorithm  to  yield  useful 
results  is  dependent  on  the  acceptable  error  and  the  signal-to-noise 
ratio.  Since  the  signal-to-noise  ratio  is  a  function  of  the  reflection 
coefficient  of  the  panel  the  error  associated  with  a  particular  time 
window  will  vary  from  one  panel  to  another.  In  order  to  obtain  some  idea 
of  the  minimum  time  window,  the  steel  panel  described  earlier  was 
investigated. 

Three  waveforms  were  obtained  for  the  reflection  of  a  step  sinusoid 
from  the  steel  panel  at  frequencies  of  2,  2.5,  and  3  kHz.  The  signal- 
to-noise  ratio  in  each  was  approximately  40  dB.  Each  waveform  was 
analyzed  with  a  15-pole  expansion  that  included  two  steady-state  driving 
frequency  poles  entered  a  priori.  The  time  window  was  varied  from  250 
to  130  ps  in  10-ps  steps. 


Fig.  (9)  -  Amplitude  of  the  steady-state  driving 

frequency  poles  as  a  function  of  time  window 
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In  Fig.  (9)  the  amplitude  of  the  driving  frequency  pole  has  been 
plotted  against  the  length  of  the  time  window  for  each  of  the  waveforms. 
Since  the  actual  amplitudes  are  unknown,  the  average  values  of  the  flat 
portions  of  the  curves  were  used  as  the  correct  amplitudes .  Arbitrarily 
choosing  an  allowable  error  of  0.25  dB  from  the  average  value  as  a  measure 
of  the  accuracy  of  the  algorithm  produced  minimum  time  windows  of  140  ys 
for  the  3-kHz  waveform,  170  ys  for  the  2.5-kHz  waveform,  and  230  ys  for 
the  2-kHz  waveform.  These  time  windows  correspond  to  0.42,  0.425,  and 
0.46  wavelengths,  respectively,  for  2.5  and  2  kHz.  Since  these  values 
are  in  good  agreement,  a  general  rule  of  approximately  half  a  wavelength 
as  the  minimum  time  window  has  been  used  for  the  data  obtained  in  this 
report . 

EXPERIMENTAL  PROCEDURE 

There  are  two  procedures  for  analyzing  echo-reduction  measurements 
with  the  modified  Prony  method.  In  the  first,  referred  to  as  the  two- 
window  method,  the  projector  is  positioned  170  cm  and  the  probe  15  cm 
in  front  of  the  panel.  A  USRD  type  F36  standard  transducer  is  used  as 
the  projector  while  the  probe  is  a  USRD  type  H52  standard  hydrophone. 

The  projector  is  driven  by  a  step  sinusoidal  signal  of  1-ms  duration 
that  produces  a  200-ys  segment  of  incident  signal,  at  the  probe, 
followed  by  200  ys  of  incident  plus  reflected  signal  before  the  arrival 
of  the  diffracted  signal  from  the  panel  edges.  This  allows  equal  periods 
of  the  incident  and  reflected  signals  to  be  observed. 

The  waveform  at  the  probe  is  sampled  at  1  MHz,  and  approximately 
50  to  100  waveforms  are  averaged  to  reduce  the  incoherent  noise  level. 

The  waveform  is  then  divided  into  two  time  windows — one  containing 
only  the  incident  signal,  and  the  second  containing  the  incident  plus 
reflected  portions  of  the  signal.  Figure  (10)  illustrates  the  waveform 
and  the  two  time  windows  used  in  analyzing  the  waveform.  Both  time 
windows  are  analyzed  by  the  modified  Prony  method  to  find  the  amplitudes 
of  the  steady-state  driving  frequency  poles.  A  15-pole  expansion  is  used 
with  three  poles  entered  a  priori.  Two  of  the  a  priori  poles  are  the 
complex  conjugate  pair  representing  the  driving  frequency,  and  the  third 
a  priori  pole  is  a  real  pole  associated  with  a  high-pass  RC  filter  on  the 
input  side. 
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Fig.  (.10)  -  Illustration  of  two-window  method.  Two 

time  windows  are  obtained  from  the  waveform; 
the  first  concaining  only  the  incident  signal 
and  the  second  containing  the  incident  plus 
reflected  signal. 

The  complex  amplitude  of  the  incident  portion  is  then  phase  shifted 
by  an  amount  equal  to  the  time  separation  of  the  two  time  windows  and 
substracted  from  the  complex  amplitude  of  the  incident  plus  reflected 

This  yields  the  complex  amplitude  of  the  steady— state  driving 
frequency  pole  for  the  reflected  signal.  The  echo  reduction  is  then 
calculated  as 


Echo  reduction  =  20  log  /A 


(35) 


where  Aj  and  Ar  are  respectively  the  moduli  of  the  amplitudes  of  the 
incident  and  reflected  signals. 

While  the  two-window  method  yields  good  results,  it  is  not  as 
accurate  as  the  second  method,  to  be  discussed  below,  owing  to  phase 
errors.  The  algorithm  does  a  much  better  job  of  finding  the  correct 
modulus  of  the  amplitude  than  it  does  in  finding  the  correct  phase. 
This  phase  error  introduces  an  error  into  the  echo-reduction  calcu¬ 
lation  when  the  incident  amplitude  is  phase  shifted  and  substracted 
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from  the  incident  plus  reflected  amplitude.  Not  all  of  the  incident 
signal  is  cancelled,  and  the  amplitude  obtained  does  not  represent  the 
reflected  signal  only.  The  magnitude  of  the  error  will  depend  on  the 
phase  error  and  the  relative  phase  of  the  direct  and  reflected  signals. 
If  the  length  of  the  time  window  is  equivalent  to  at  least  one  period 
of  the  driving  frequency,  frequencies  above  5  kHz  for  a  200-ys  time 
window,  the  phase  error  is  negligible  and  the  two-window  method  is 
sufficient.  However,  as  the  frequency  is  reduced,  the  phase  error 
increases  and  an  alternate  method  must  be  vised. 

The  second  method,  referred  to  as  the  difference  method,  eliminates 
the  effect  of  the  phase  error  in  the  algorithm  by  directly  subtracting 
out  the  incident  signal.  This  method  was  basically  described  in 
connection  with  data  acquisition  for  the  section  on  the  data  increment. 
It  consists  of  performing  two  separate  measurements — one  with  and  one 
without  the  acoustic  panel  in  position.  The  two  recorded  waveforms 
are  then  directly  subtracted  to  yield  the  reflected  signal  from  the 
panel  as  illustrated  in  Fig.  (11) .  Then  the  modified  Prony  method  is 
used  to  analyze  the  incident  and  reflected  waveforms  separately  to 
obtain  the  amplitudes  of  the  driving  frequency  poles. 
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Fig.  (11)  -  Sequence  illustrating  difference  method.  The 

top  waveform  was  obtained  with  panel  in  position 
while  the  middle  waveform  was  obtained  without  the 
panel.  The  lower  waveform  was  obtained  by  direct 
substraction  of  the  width  and  without  panel  waveforr 


I 


The  difference  method  has  the  additional  advantage  of  allowing 
longer  portions  of  the  reflected  signal  to  be  observed.  In  the  two- 
window  method  the  optimum  position  of  the  probe  is  15  cm  in  front  of 
the  panel  since  this  allows  equal  segments  of  the  incident  and  reflected 
signals  to  be  observed.  However,  in  the  difference  method  the  measure¬ 
ment  performed  without  the  panel  produces  the  waveform  for  the  incident 
signal.  This  allows  the  panel  to  be  positioned  close  to  the  probe  for 
the  second  measurement.  With  the  probe  5  cm  from  the  panel,  250  us  of 
reflected  signal  can  be  observed  prior  to  the  arrival  of  the  diffracted 
signal  from  the  panel  edges.  However,  care  must  be  taken  to  insure  that 
the  two  measurements  are  identical.  In  addition  a  least-square  sub¬ 
traction  should  be  used  in  obtaining  the  reflected  signal  to  compensate 
for  any  gain  and  phase  changes  that  may  occur  between  measurements. 

The  disadvantages  of  the  difference  method  are  the  additional  time 
required  for  separate  measurements  and  an  inherent  phase  error  due  to 
digitizing  the  waveform.  The  time  factor  essentially  doubles  the  time 
required  to  perform  the  measurements  while  the  phase  error  becomes  a 
problem  only  at  high  frequencies  where  the  two-window  method  is  accurate. 
This  results  in  an  obvious  choice  of  using  the  two-window  method,  except 
at  low  frequencies  (below  5  kHz)  where  the  difference  method  is  more 
accurate . 

EXPERIMENTAL  RESULTS 

In  Figs.  (12) ,  (13) ,  and  (14)  the  results  of  echo-reduction  measure¬ 
ments  of  steel  and  aluminum  panels  have  been  plotted  against  theoretical 
curves.  The  measurements  were  performed  on  0.95-cm-thick  and  76.2-cm- 
square  panels  in  the  anechoic  tank  at  USRD.  The  noise  level  in  the 
anechoic  tank  during  the  measurements  was  approximately  40  dB  below  the 
incident  signal  level. 

In  Fig.  (12)  the  measurements  of  the  steel  panel  were  processes  by 
the  difference  method.  Measurements  were  made  with  and  without  the  panel 
in  position,  and  35  waveforms  at  each  frequency  were  averaged.  The  no¬ 
panel  waveforms  were  directly  subtracted  from  the  waveforms  with  the 
panel  in  position  to  obtain  the  waveform  of  the  reflected  signal.  Each 
waveform  was  then  processed  with  a  15-pole  expansion  that  included  the 
two  driving  frequency  poles  a  priori.  The  results  deviate  from  the 
theoretical  curve  by  a  few  tenths  of  a  dB. 
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ECHO  REDUCTION  (dB) 


In  Figs.  (13)  and  (14)  the  results  of  measurements  on  an  aluminum 
panel  have  been  plotted  against  theoretical  curves.  The  aluminum 
plate  was  chosen  since  it  has  a  larger  echo  reduction  than  the  steel 
panel  and  was  a  better  test  of  the  method.  The  results  in  Fig.  (13) 
were  obtained  in  the  same  manner  as  the  data  for  the  steel  panel 
except  three  a  priori  poles  were  used.  The  third  a  priori  pole  was 
associated  with  an  RC  filter  on  the  input  side  of  the  electronics. 

The  experimental  results  deviate  from  the  theoretical  curve  by 
approximately  0.25  dB. 

In  Fig.  (14)  the  measurements  were  processed  by  the  two-window 
method.  The  probe  was  positioned  15  cm  from  the  panel  to  provide 
equal  segments  of  the  incident  and  reflected  signals  and  50  separate 
measurements  were  averaged  to  obtain  the  waveforms.  Each  time  window 
was  processed  with  a  15-pole  ejqpansion  that  included  three  a  priori 
poles.  The  results  deviate  from  the  theoretical  curve  by  an  average 
of  0.35  dE.  However,  the  3-kHz  measurement  deviates  by  1.05  dB,  a 
result  explained  by  the  previously  described  phase  error  associated 
with  the  two-window  method. 

In  addition  to  the  data  presented  here,  measurements  on  a  1.27- 
cm-thick  steel  plate  have  been  performed  at  2  and  2.5  kHz  with  the 
difference  method.  These  measurements  deviated  from  the  theoretical 
values  by  approximately  0.5  dB  and  indicate  that  the  method  is  capable 
of  performing  accurate  measurements  down  to  2  kHz. 
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Fig.  13  -  Results  of  Prony 
measurements  of  a  0.95-cm- 
thick  aluminum  panel  processed 
by  the  difference  method. 


Fig.  14  -  Results  of  Prony 
measurements  of  a  0.95-cm- 
thick  aluminum  panel  processed 
by  a  two-window  method. 
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CONCLUSION 


The  experimental  results  indicate  that  the  modified  Prony  method 
is  capable  of  making  echo-reduction  measurements  down  to  2  kHz  on 
simple  homogeneous  panels  with  an  error  no  greater  than  0.5  dB.  There 
have  been  no  measurements,  as  yet,  on  high  Q  or  lossy  panels.  However, 
these  panels  should  not  present  an  obstacle  as  long  as  the  required 
number  of  terms  in  the  expansion  does  not  become  too  large. 
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APPENDIX  A 


OVERVIEW 

The  following  modified  Prony  program  listing  has  been  written  in 
FORTRAN  4+  and  is  compatible  with  the  Digital  Electronics  Corporation 
PDP  11/45  computer  with  the  system  RSX-11D.  The  program  is  designed 
to  use  data  files,  with  a  maximum  of  1024  data  points,  that  have  the 
short  IAG  header  format.  Current  dimension  statements  have  limited 
the  program  to  a  maximum  order  of  expansion  of  15  and  a  maximum  of 
5  a  priori  poles. 
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APPENDIX  B 


LIST  OP  VARIABLES  IN  PRONY 


J 

- 

number  of  data  points  in  data  file 

ISP 

- 

initial  start  point  in  data  file 

NPTS 

- 

number  of  points  in  time  window 

IBD 

- 

data  increment 

NR 

- 

order  of  expansion 

IA 

- 

number  of  a  priori  poles 

DX 

- 

data  file  time  increment 

DEV 

- 

mean  square  deviation 

ROOTS (I) 

- 

array  of  a  priori  poles  (s  plane) 

ROOTZ(I) 

- 

array  of  z  plane  poles 

DATA (I) 

- 

array  of  data  points 

COE (I) 

- 

coefficients  of  Prony  polynomial 

ACDEF(I) 

- 

array  of  amplitudes 

COEFB(I) 

- 

coefficients  of  a  priori  polynomial 

E(I,K) 

- 

matrix  of  Prony  difference  equations 

F(I) 

vector  associated  with  Prony  difference 

equations 

A(I,K) 

- 

matrix  of  equations  for  amplitudes 

APPENDIX  C 


INPUT  EXAMPLE 

RUN  PHONY  ) 

ENTER  FILE  SPECS ♦ 

ERON.DHTillS 

ENTER  INITIAL  START  POINT ♦ 25? 3 
ENTER  NUMBER  OF  POINTS  IN  TIME  UINDQU.200  ^ 
ENTER  BASIC  DATA  INCREMENT . (ODD  INTEGER)  11  ) 
ENTER  ORDER  OF  EXPANSION.  15.^ 

DO  YOU  WISH  TO  ENTER  APRIORI  ROOTS?  (Y/N>  1 } 
ENTER  t  OF  APRIORI  POLES’?^ 

ENTER  POLE  VALUES  AS 

REAL  >  I M A G I N A R Y 


1  1^5  55.2  ) 

2  <L^i-l8STVt55^?  ) 

FINISHED 


Note:  I) ,  All  underlined  portions  are  user  supplied 
2),  )  indicates  ’RETURN’ 


29 


APPENDIX  D 


OUTPUT  EXAMPLE 

DATE®  03- MAR- 80 

TIME®  10 s 30: 03 

INPUT  DATA  FILE  INFORMATION 

PRON.DHT; 11 

FILE  HEADER  INFORMATION 

o  1024  i  e.oeeoo  e.ieeeeE-03  c  2 

3  KHZ 

INITIAL  START  POINT®  252 

TIME  WINDOW*  0.20000E-03 
BASIC  DATA  INCREMENT®  11 

ORDER  OF  EXPANSION®  15 
APRIORI  POLES 
*  REAL.  1MAC 


1  0.OOOOOE+OO  0. 1CC50E+03 

2  0.OOCO3E+OQ-O. 1C350E+03 

Z- PLANE  POLES.  RESIDUES. 


1 

0.999G2E+00 

0. 1884GE-01 

1 

0.00235E-01 

-0. 10084E+00 

2 

0.999C2E+QQ 

-0. 18C48E-0 I 

o 

0.60255E-01 

0. 1G0C4E+C0 

3 

O.99G93E+O0 

0. 62062E-0 I 

3 

-0.3704  IE-03 

0 . 50 109E-03 

4 

0.99C93E+00 

-0 . 62062E-01 

4 

-0.37040E-03 

-0 . 50 1CCE-03 

5 

0.  10032E+0 1 

0.  101 19E+00 

5 

-0 . 3703QE-03 

0.42420E-03 

6 

0.  1C032E+01 

-0. 101 19E+00 

6 

-0. 37030E-03 

-0. 42420E-03 

7 

0 . 97 I39E+00 

0.  14830E+00 

7 

-0.3 1909E-02 

-0.21330E-03 

0 

0.97 139E+00 

-0. 14GS0E+00 

8 

-0 . 3  19G9E-02 

0 . 2 1537E-03 

0 

0. 90003 E+00 

0.  13233E+00 

9 

-0 . 47842E-03 

-0. 60498E-03 

10 

0.98065E+00 

-0. 18235E+00 

10 

-0. 47042E-03 

0. 60498E-O3 

1 1 

0. 96659E+Q0 

0 . 235G0E+00 

11 

0. 16059E-03 

0.74164E-04 

12 

0.96639E+00 

-0.25580E+00 

12 

0. 16059E-03 

-0 . 74 165E-04 

13 

0.9.T033F.+00 

0.20486E+00 

23 

0. 15642E-02 

-0. 1331CL-02 

14 

0.930330+09 

-0.2040&F.+00 

14 

0. 130  12E-02 

0.  13310E-02 

MEAN  SQUARE  DEVIATION*  0.30100E-C6 


o  onn 


APPENDIX  E 


PRONY  LISTING 


FRONY  MAIN  PROGRAM 

COMPLEX  ROQTSUO) > ROOTZ < 15 ) , ACOEF < 15) »E(15»15)»F(15) 
DIMENSION  COEFB< 16) , A ( 15 » 15 ) » R< 15 ) » DATA ( 1024 ) , COE ( 16 ) 
BYTE  TIM(S) »  DAT ( 9 ) 

C 

CALL  DATE < DAT) 

WRITE(2»500)DAT 
CALL  TIME ( TIM ) 

WRITE  <2»510)TIM 
CALL  RFIL.ES  <  DATA  » J » DX ) 

C 

C  OBTAIN  INITIAL  START  POINT 

C 

1  WRITE (6, 520) 

READ (5, 530) ISP 
WRITE<2»535) ISP 
C 

C  CHECK  THAT  ISP  IS  GREATER  THAN  ZERO 

C 

IF< ISP-1 )5» 10. 10 
5  WRITE (6 » 540) 

GOTO  1 
C 

C  OBTAIN  NUMBER  OF  POINTS  IN  DATA  WINDOW 

C 

10  WRITE  <  6 » 550 ) 

READ ( 5 >  560 ) NPTS 
TW~DX*NPTS 
C 

C  CHECK  THAT  WINDOW  DOES  NOT  EXCEED  DATA  RANGE 

C 

IF  < ISP+NPTS- J-l ) 20  *  20 » 15 
15  WRITE ( 6 1 570 ) 

GOTO  1 

20  WRITE (2* 580) TW 

C 

C  OBTAIN  BASIC  DATA  INCREMENT 

C 

25  WRITE  <  6 r 590 ) 

READ <5» 600) IBD 


|  PHSC*0U» 


c 

C  CHECK  THAT  IBD  IS  AN  ODD  INTEGER 

C 

I J=IBD/2 

IF  <  2#I J-IBO ) 35, 30 »  30 
30  WRITE  <6,610) 

GOTO  25 
C 

C  OBTAIN  ORDER  OF  EXPANSION 

C 

35  WRITE <6 , 620) 

READ <  5 , 630 ) NR 
WRITE ( 2  f  640) NR 
C 

C  ARE  THERE  ANY  APRIORI  ROOTS? 

C 

WRITE(6»650) 

READ <5, 660) 12 
IFdZ-'Y'  )40,50,40 
C 

C  NO  APRIORI  ROOTS 

C 

40  I A=0 

COEFB( 1 )=1 .0 
COEFB( 2 ) =0 . 0 
COEFB ( 3 ) =0 . 0 
C 

C  CHECK  THAT  THERE  ARE  SUFFICIENT  DATA  POINTS 

C 

IF ( NPTS-NR# ( IBD+1 > >45,70,70 
45  WRITE ( 6 » 670 ) 

GO  TO  10 

C 

C  CALL  PRIORI  FOR  APRIORI  ROOTS 

C 

50  CALL  PRIORI < ROOTS » ROOTZ  »  COEFB » I A  f  DX  r IBD ) 

C 

C  CHECK  THAT  IA  WAS  NOT  SET  EQUAL  TO  ZERO  IN  PRIORI 

C 

IF ( IA ) 55  >40*55 
C 

C  CHECK  THAT  THERE  ARE  SUFFICIENT  DATA  POINTS  FOR  THE 

C  CASE  WITH  APRIORI  ROOTS. 

C 

55  IF < NPTS-NR* < IBD+l )+I A >60,65,65 

60  WRITE  <  6 , 670 ) 

GOTO  10 
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CHECK  FOR  CASE  WHERE  NUMBER  OF  APRIORI  ROOTS= 
ORDER  OF  EXPANSION 


C 
C 
C 
C 

65  IF(NR-IA)66» 110,70 

66  NR=IA 
GOTO  110 

70  CALL  MATRIX  <  A » R  *  NR » IA  >  NPTS  r COEFB  r DAT  A  > ISP » IBD  > IL ) 

C 

C  CHECK  RETURN  FROM  MATRIX 

C 

IF(IL-10)75>10»75 

C 

C  CHECK  FOR  THE  CASE  WHERE  THERE  IS  ONLY  ONE  NON-APRIORI 

C  ROOT 

C 

75  IF(NR-IA-1)65»80>85 

C 

C  OBTAIN  SINGLE  DESIRED  ROOT 

C 

80  R00T2<IA+1)=CMPLX(-R(1)/A(1,1> ,0.0) 

GOTO  110 

85  CALL  SOLVER  <  A  *  R  >  COE  *  NR  r I A  > I B ) 

C 

C  CHECK  EXIT  FROM  SOLVER 

C 

IF ( IB ) 95  f  90 »  95 
90  NR=NR-1 

WRITE ( 2  j  680 ) NR 
IF ( NR-IA ) 65 » 1 10 »  70 

95  CALL  F'RQD  <  COE » NR  >  I ER >  ROOTZ » I A ) 

C 

C  CHECK  RETURN  FROM  PROD 

C 

IF  ( IER )  1 1 0 » .1 10  » 105 
105  NR=NR-1 

GO  TO  70 

.1.10  IF  ( I BD- 1 )  25 » 120  >  1 15 

C 

C  CONVERT  Z-PLANE  ROOTS  TO  IBP=1 

C 

.1.15  CALL  ROOTC  (  ROOTZ  » IBD  r  NR  ) 

120  WRITE  <  2 »  700 ) 

DO  125  1  =  1 » NR 

WRITE (2»710)If ROOTZ (  I ) 

125  CONTINUE 

C 


C  OBTAIN  ACOEF 
C 

130  CALL  RES  I  DU  ( DATA  >  E>  F »  ROOTZ  > NR  t  NPTS  > ISP ) 

CALL  SOLVE (E>  F > ACOEF >NR> IB) 

IF(IB)140>135>140 
135  CALL  ROOTE  <  ROOTZ  > NR ) 

GOTO  130 

140  WRITE ( 2  >  720 ) 

DO  150  1=1 >NR 

WRITE ( 2  >  710 )  I >  ACOEF ( I ) 

150  CONJINUE 

CALL  MSDEV < ISP > DATA > DEV > NPTS > ROOTZ > ACOEF > NR ) 

WRITE ( 2 »  730 ) DEV 
170  WRITE ( 6 » 750 ) 

C 

C  FORMAT  STATEMENTS 
C 

500  FORMAT  ( / > 16X  > '  DATE=  '.9A1) 

510  FORMAT  < / » 16X > '  TIME=  ' >  8A1 ) 

520  FORMAT (/ t ' CENTER  INITIAL  START  POINT.') 

530  FORMAT  < 14) 

535  FORMAT  ( / > 10X  > '  INITIAL  START  POINT=  '.14) 

540  FORMAT ( /  >  '  INITIAL  START  POINT  MUST  BE  GREATER  THAN  ZERO. 
550  FORMAT ( /  >  ' $ENTER  NUMBER  OF  POINTS  IN  TIME  WINDOW.') 

560  FORMAT (14) 

570  FORMAT (/>'  DATA  WINDOW  EXCEEDS  DATA  FILE.') 

580  FORMAT ( / > 16X  > '  TIME  WINDOW*  '»E12.5> 

590  FORMAT (/> '$ENTER  BASIC  DATA  INCREMENT .( ODD  INTEGER)  ') 
600  FORMAT (13) 

610  FORMAT ( /  >  '  BASIC  DATA  INCREMENT  MUST  BE  ODD  INTEGER.') 

620  FORMAT ( /  >  ' CENTER  ORDER  OF  EXPANSION,') 

630  FORMAT ( 13 ) 

640  FORMAT  <  7  > 16X> '  ORDER  OF  EXPANSION*  '>13) 

650  FORMAT ( /  >  ' $D0  YOU  WISH  TO  ENTER  APRIORI  ROOTS?  (Y/N)  ') 
660  FORMAT ( 1A1 ) 

670  FORMAT ( /  >  '  INSUFFICIENT  NUMBER  OF DATA  POINTS,') 

680  FORMAT ( / > 1 6X  > '  ORDER  OF  EXPANSION*  '>13) 

700  FORMAT (/> 10X> '  Z-PLAME  POLES. '»/> 

710  FORMAT (/>6X> 12 > 6X > El  2 . 5 > 5X > El  2 , 5 > 

720  FORMAT (/ > 10X , '  RESIDUES .'>/ ) 

730  FORMAT (/» 16X> '  MEAN  SQUARE  DEVIATION*  '>E12.5) 

750  FORMAT ( / > '  FINISHED') 

CALL  EXIT 
END 


u 

C  PRONY  SUBROUTINE  FOR  ENTERING  APRIORI  POLES  INTO  VECTOR 
C  OF  S-PLANE  POLES. 

C  PROGRAM  ALSO  CONVERTS  S-PLANE  POLES  TO  Z-PLANE  POLES 
C  AND  COMPUTES  B  COEFFICIENTS  FOR  PRONY  DIFFERENCE  EQUATIONS. 
C  DIMENSIONING  HAS  LIMITED  SUBROUTINE  TO  5  POLES. 

C 

SUBROUTINE  PRIORI < ROOTS > ROQTZ > COEFB > I A > DX > IBP ) 

C 

C  RQOTS= VECTOR  OF  S-PLANE  POLES 
C  ROOTZ=VECTOR  OF  Z-PLANE  POLES 

C  COEFB=VECTOR  OF  B  COEFFICIENTS  ORDERED  FROM  LOU  TO  HIGH 
C  I A=NUMBER  OF  APRIORI  ROOTS 
C  DX=DAT A  FILE  TIME  INCREMENT 
C  IBD=BASIC  DATA  INCREMENT 
C 

COMPLEX  R00TS<5) »ROOTZ< 15) >BC6) fCl 
DIMENSION  COEFB  < 6 ) 

C 

C  ENTER  NUMBER  OF  ROOTS 
C 

WRITE(6f 100) 

READ<5, 110) IA 
IF ( IA ) 90 » 90  >  5 

C 

C  ENTER  ROOTS  AS  COMPLEX  NUMBERS 
C 

5  WRITE (6, 120) 

WRITEC6» 130) 

DO  10  1=1 »IA 
WRITE <6*140)1 
READ  ( 5>>  150 )  ROOTS  ( I  ) 

10  CONTINUE 

URITEC2, 160) 

WRITE ( 2  r 170 ) 

DO  15  IX=1 > I A 

WRITE (2» 180) IX » ROOTS (IX) 

15  CONTINUE 
C 

0  CONVERT  TO  Z-PLANE 
C 

C=IBD*DX 

C1=CMPLX(C»0.0) 

DO  20  IX= 1 » I A 
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i  r o  o  oocnoc.no  «  o  o  cn  o 


ROOTZ  (IX)  =CEXF‘  ( ROOTS  ( I X ) *C1 ) 

CONTINUE 

CONTINUE 

COMPUTE  B  COEFFICIENTS ,  SET  ALL  TERMS  IN  B(I)=  1 

DO  30  1=1,6 

B(I) =CMPLX  < 1 ♦ 0 , 0  *  0 ) 

CONTINUE 

CHECK  VALUE  OF  IA  AND  SET  INITIAL  VALUES 

IF ( IA-1 )90>40»45 
COEFB < 1 ) =REAL (ROOTZ(l)  ) 

COEFB  <  2 )  =  1 . 0 
GO  TO  90 

B ( 1 ) =R00TZ  < 1 ) #R00TZ  <  2 ) 

B ( 2 ) =- (  ROOTZ  < 1 ) +ROOTZ ( 2 ) ) 

IF(IA-2)90»65»  50 


ENTER  LOOP  FOR  CALCULATING  COEFFICIENTS 

DO  60  K=3 » I A 

DO  55  J=K ,2,-1 

B( J)=B< J-l) -ROOTZ (K)*B< J) 

B(1)=-B<1>*R00TZ(K) 

CONTINUE 
DO  70  1=1 » IA+1 
COEFB ( I ) =REAL ( B  <  I )  ) 

CONTINUE 

RETURN 


FORMAT  STATEMENTS 


FORMAT (/, 'CENTER  *  OF  APRIORI  POLES?') 
FORMAT (12) 

FORMAT ( / , '  ENTER  POLE  VALUES  AS',/) 
FORMAT (/,  16X,  '  REAL.,  IMAGINARY',/) 


140  FORMAT  </,'$' 12 , 3X ) 


150  FORMAT (2E12.5) 

1 60  FORMAT ( / ,27X, '  APRIORI  POLES',/) 

170  FORMAT ( 7X , '  *',15X»'  REAL,'6X,'  IMAG',/) 

ISO  FORMAT  </,7X,I2,13X,2E12«5) 

END 


PRONY  SUBROUTINE  FOR  GENERATING  THE  PRONY  DIFFERENCE 
EQUATIONS  IN  THE  FORM  A*(COE>=R,  WHEN  APRIORI  POLES 
( R00T2 )  ARE  GIVEN. 'A'  IS  GENERATED  BY  A  VIRTUAL  MATRIX 
PREMULTIPLIED  BY  IT'S  TRANSPOSE. 

SUBROUTINE  MATRIX. ( A > R . NR > I A > NPTS. CQEFB > DATA t ISP > IBP t IL ) 

A=MATRIX  CONTAINING  DIFFERENCE  EQUATIONS 
R-COLUMN  VECTOR  WHICH  ARISES  DUE  TO  THE  CONSTRAINT 
THAT  THE  HIGHEST  BETA  COEFFICIENT  EQUAL  ONE 
NR=NUMBER  OF  POLES  IN  PRONY  EXPANSION 
IA=NUMBER  OF  APRIORI  POLES 
NPTS=NUMBER  OF  POINTS  IN  DATA 
COEFB= COEFFICIENTS  FROM  APRIORI  POLES 
DATA=DATA  FILE 

ISP= INITIAL  START  POINT  IN  DATA 
IBD=BASIC  DATA  INCREMENT 
IL=RETURN  CODE 

DIMENSION  A(15j15)»R<15) t COEFB (6) t DATA ( 1024 ) 

DEFINE  VARIABLE  RANGE 

IC=NR-IA 

IR=NPTS-NR*IBD 

IS=ISP-1 

CHECK  FOR  EXACTLY  SOLVED  CASE 

IN=IR-IC 

IF< IN) 10.20.50 

IL-10 

WRITE ( 6 t 500 ) 

RETURN 

EXACTLY  SOLVED  CASE 

DO  30  1=1 rIC 
DO  30  J=I.IC 
A(I. J)=0.0 
DO  30  I K  =  1  f  I A  + 1 


% 


A(I  t  J)=COEFB(IK)*DATA(  I  +  IS+dK-1  >*IBD+(  J-1)*IBD)+A<  I ,  J) 

30  CONTINUE 

DO  40  1  =  1  dC 

R ( I ) =0 . 0 

DO  40  IK=1 f IA+1 

R(I)=R( I )-C0EFB< IK)*DATA< I  +  IS+ « IK-1 ) *IBD+IC*IBD ) 

40  CONTINUE 
GO  TO  SO 
C 

C  LEAST  SQUARE  TYPE  SOLUTION 
C 

50  DO  70  1=1 fIC 
DO  70  J=I»IC 

A(Ii J)-0.0 
R  ( I ) =0 « 0 
DO  70  K'=1 » IR 
B1 =0 . 0 
B2=0 . 0 
R1=0 « 0 

DO  60  IK=1 » IA+1 

Bl=COEFB<IK>*DATA<K+IS+(IK-l>#IBD+< J-1)*IBD)+B1 
B2=CQEFB< IK)*DATA(K+IS+< IK-1 ) *IBD+ < 1-1 > *IBD ) +B2 
Rl=COEFB(IK)*DATA<K+IS+< IK-1 ) *IBD+IC*IBD > +R1 
60  CONTINUE 

A(Ii J)=B1*B2+A(I, J) 

R  < I ) =R < I ) -ftl #B2 
70  CONTINUE 
80  IF(IC-l)110fll0>90 

90  DO  100  1  =  2  iIC 
DO  100  J=1»I-1 
A(I»J)=ACJ»I) 

100  CONTINUE 
110  IL=36 
RETURN 

500  FORMAT ( / t '  INSUFFICIENT  NUMBER  OF  DATA  POINTS»IN  (MATRIX).') 
END 
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PRONY  SUBROUTINE  FOR  SOLVING  THE  LEAST  SQUARE  EQUATIONS 
GENERATED  IN  MATRIX  TO  FIND  THE  COEFFICIENTS  OF  THE 
PRONY  POLYNOMIAL 


SUBROUTINE  SOLVER ( A » R f COE > NR > I A » IB ) 


DIMENSION  A(15»15)»R(15) fC0E<16)fX(15 

IB=1 
N=NR-IA 
DO  10  1=1 *N 
IK'TA  <  I )  =1 
CONTINUE 
K=1 

CHECK  LEADING  TERM 

I F ( A ( K »  K ) ) 30 »  20  ?  30 

CALL  INTERD (A>R»IKTA>K»NrICl) 

IF(IC1)30,25,30 

IB=0 

RETURN 

CONTINUE 

DIVIDE  ROUS  BY  LEADING  TERM 


C1=A (K *  K ) 

R(K)=R(K)/C1 
DO  40  J=K  f N 
A(K» J)=A(Kr J)/C1 
CONTINUE 

SUBTRACT  K  ROW  FROM  ALL  ROWS  BELOW 

DO  50  I =K+ 1 1 N 

R< I )=R( I>-R(K)*A(I rK) 

Cl =A ( I »  K ) 

DO  50  J=K » N 

A<I,  J)=A<I> J)-A(K> J)*C1 
CONTINUE 


)  t IKTA(IS) 


K=K+1 

IF ( N-N ) 15 » 60 »  60 
60  X(N)=R<N)/A<N*N) 

DO  70  I-N-l » 1 » -1 
X  < I ) =R  < I ) 

DO  70  J=N» 1+1 » -1 
X<I)=X<I)-X<J>*A<I»J> 
70  CONTINUE 

DO  80  1=1 fN 
J=IKTA< I ) 

COE< J)=X( I ) 

80  CONTINUE 

COE  (  N+l )  =  1 . 0 

RETURN 

END 
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no  non 


C 

C 

SUBROUTINE  INTERCHANGES  ROWS  AND  COLUMNS  OF  MATRIX  A  AND 
VECTOR  R  WHILE  KEEPING  TRACK  OF  CHANGES  IN  VECTOR  IKTA . 

SUBROUTINE  INTERD< At R» IKTAt K»N> IC1 ) 


DIMENSION  A<15»15)»R<i5)»XlT(15)» X2T <15)»IKTA(16) 
C 
C 

IC1=0 
IR1  =K 
X1=0.0 
DO  10  I=KfN 

IF ( ABS (XI )-ABS(A<K»I) >)5f10f10 
5  IC1=I 

X1=A (K» I ) 

.10  CONTINUE 

IF ( ICl ) 15f 15  f  20 
15  RETURN 

20  I«IR1 

DO  30  J=1 f  N 
X1T< J)=A< I f J) 

30  CONTINUE 

I  =  IC1 

DO  40  J=1fN 
X2T( J)=A< If J) 

40  CONTINUE 

I  =  I R 1 

DO  50  J= 1 1  N 
A  (  I t J ) =X2T ( J ) 

50  CONTINUE 

I  =  ICl 

DO  60  J= 1 «  N 
A  <  I  f  J  )  =  X 1  T  \  J  ) 

60  CONTINUE 

C 

C  INTERCHANGE  COLUMNS 

C 

J=[R1 

DO  70  1=1 fN 
X 1 T  < I ) =A (  I  f  J ) 
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70 

SO 

90 

110 


CONTINUE 

J=IC1 

DO  SO  1  =  1, N 
X2T<I)=A<I, J) 
CONTINUE 
J=IR1 

DO  90  1=1, N 
A ( I , J) =X2T ( I ) 
CONTINUE 
J=IC1 

DO  110  1=1, N 
A(I, J)=X1T(I) 
CONTINUE 
R1T=RUR1) 
R2T=R< IC1 ) 

R ( IR1 ) =R2T 

R<IC1)=R1T 

I=INTA ( IR1 ) 

J*INTA(IC1) 

IKTA(IR1)=j 

IKTA<IC1)=I 

RETURN 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


A 


8 

10 

20 


C 

c 

c 

c 

c 


10 
1 5 


PRONY  SUBROUTINE  FOR  CONVERTING  THE  ROOTZ  FOUND  WITH 
THE  BASIC  DATA  INCREMENT  NOT  EQUAL  TO  ONE  (ROOTZ )**IBD 
TO  THE  ROOTZ  WITH  BASIC  DATA  INCREMENT  OF  ONE. 

SUBROUTINE  ROQTC ( ROOTZ > IBP f NR ) 

ROOTZ=CONTAINS  THE  ROOTZ** IBD  ON  RETURN  CONTAINS  ROOTZ 
IBD-BASIC  DATA  INCREMENT 

COMPLEX  ROOTZ (15) 


DO  20  1=1 » NR 
A=AIMAG ( ROOTZ ( I )  ) 

B=REAL ( ROOTZ ( I ) ) 

J=1.0E+3*A 
K=ININT( 1 . 0E3*B ) 

IF(J)10*5*10 

IF(K)6»8»10 

B=AB3 ( REAL (ROOTZ(I))  > 

B=B** ( 1/FLOAT ( IBD ) ) 

ROOTZ ( I ) =CMPLX ( -B » 0 , 0 ) 

GOTO  20 

ROOTZ ( I ) =CEXP ( ( CLOG ( ROOTZ ( I ) ) ) /IBD ) 

CONTINUE 

RETURN 

END 

PRONY  SUBROUTINE  FOR  DELETING  ROOTS  WHEN  RESIDU  FAILS 

SUBROUTINE  ROOTE ( ROOTZ »  NR  > 

COMPLEX  ROOTZ (15) 


WRITE (6 v 100) 

DO  10  1=1 t NR 

WRITE ( 6  f 1 10 ) I f  ROOTZ ( I ) 

CONTINUE 

WRITE ( 6  r 1 20 ) 

READ ( 5 1 130)  IX 


45 


DO  20  I=IX>NR-1 
ROOTZ  < I ) =ROOTZ ( 1  +  1 > 

20  CONTINUE 

ROOTZ  <  NR ) =CMPLX ( 0 . 0 »  0 . 0 ) 

NR=NR-1 
WRITE (6r 140) 

READ (5* 150) I J 
IF(IJ-'r')30»15r30 
30  RETURN 

C 

C  FORMAT  STATEMENTS 

C 

100  FORMAT</» ' 1Z-PLANE  POLES') 

110  FORMAT (3X»I2»3X>E12.5f4X»E12.5) 

120  FORMAT (/ r '$UHICH  POLE  IS  TO  BE  DELETED?') 

130  FORMAT (12) 

140  FORMAT (/, ' ^DELETE  ANOTHER?') 

150  FORMAT ( 1 A1 ) 

END 
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^ mm 


n  n  n  o  o 


PRONY  SUBROUTINE  WHICH  LOADS  THE  MATRIX  E  AND  VECTOR  F 
WITH  THE  LEAST  SQUARE  EQUATIONS  FOR  CALCULATING  THE 
RESIDUES  ASSOCIATED  WITH  THE  POLES. 


SUBROUTINE  RESIDU< DATA . E  .  F . RC0T2 .  NR . NPTS .ISP) 
C 

C 

COMPLEX  ROOTZ (15) »A.B.EC15»15).F<15) 

DIMENSION  DATA (1024) 

C 

C 

A=CMPLX( 1 .0.0.0) 

DO  10  1=1. NR 
DO  10  J=I.NR 

B=  <  CON JG ( ROOTZ ( I ) ) ) *ROOTZ ( J ) 

IF(AIMAG(B) )5.1»5 

1  IF ( REAL (B)-1.0)5»2>5 

2  Ed.  J)=CMPLX(  FLOAT  (NPTS)  .0.0) 

GO  TO  10 

S  E(I. J)=(A-B**(NPTS) >/(A-B) 

10  CONTINUE 

C 

C 

DO  20  K— 1 . NR 

F(K)=CMPLX( DATA (ISP) .0.0) 

DO  20  1=1. ( NPTS- 1) 

F ( K ) =DATA ( ISP+I ) * <  CON JG ( ROOTZ ( K >  >  ** I >  +F (K ) 

20  CONTINUE 

DO  30  1=2. NR 

DO  30  J=1.I--1 

Ed.  J)=CONJG(E(J.I)  ) 

30  CONTINUE 

RETURN 
END 
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PRQNY  SUBROUTINE  WHICH  SOLVES  THE  LEAST  SQUARE 
EQUATIONS  GENERATED  IN  RESIDU  TO  FIND  THE  RESIDUES. 
SUBROUTINE  SOLVES  SIMULTANEOUS  EQUATIONS  WITH 
COMPLEX  COEFFICIENTS. 

SUBROUTINE  SOLVE (E,F,ACOEF. NR.  TFt) 


COMPLEX  E<15»15)iF(15)» ACOEF <15),X(15)»C1»C2 
DIMENSION  IKT ( 15 ) 


FILL  ARRAY  TO  KEEP  TRACK  OF  ROW  AND  COLUMN 
INTERCHANGES 

DO  10  1=1, NR 
INT< I )=I 
CONTINUE 
K=1 

CHECK  LEADING  TERM 


C3=REAL ( E  <  K , K ) ) 

IF ( C3 ) 34 , 20 » 34 

INTERCHANGE  ROWS  AND  COLUMNS 

CALL  INTERC (E, F , IKT, K, NR, IC1 ) 

IF(IC1 >34,30,34 

IB=0 

RETURN 

DIVIDE  ROWS  BY  LEADING  TERM 
N=K 

C2=£ ( N , N ) 

F(N)=F(N)/C2 
DO  40  J=K , NR 
E(N«J)=E(N,J) /C2 
CONTINUE 


cn  n  o  ri 
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IF(NR-N)50f50f45 

N=N+1 

GO  TO  35 


SUBTRACT  K  ROW  FROM  ALL  ROUS  BELOW 


0  DO  60  I=K+1 t NR 

F<I)=F<I)-F(K)*E(I»K) 
C2=E< I fK> 

DO  60  J=K>NR 
E(If J)=E<I» J)-E<Kt J)*C2 
60  CONTINUE 

K=K+1 

IF  <  K-NR) 15  >  70 » 70 
70  X ( NR ) =F  <  NR ) /E  <  NR  *  NR ) 

DO  80  I=NR-1 » 1 » -1 
X(I)=F(I) 

DO  30  J=NR » 1+1 » -1 
X<I)=X<I)-X<J)*E<I> J> 

80  CONTINUE 

DO  90  1=1 f NR 
J=IKT ( I ) 

ACOEF ( J ) =X ( I ) 

90  CONTINUE 

RETURN 
END 
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no  no  noooo 


SUBROUTINE  INTERCHANGES  ROWS  AND  COLUMNS  OF  MATRIX  E 
AND  VECTOR  F  WHILE  KEEPING  TRACK  OF  CHANGES  IN  VECTOR 
IKT. 

SUBROUTINE  INTERC (E > F > IKT > K > NR > IC1 ) 


COMPLEX  E( 15,15) »F(15) >X1T<15) »X2T<15) »R1T»R2T 
DIMENSION  IKT ( 1 5) 


IC1=0 

IR1=K 

X1=0.0 

DO  10  I=K » NR 
X2=REAL ( E ( K » I ) ) 

IF ( ABS ( XI ) -ABS  <  X2 ) )5»10»10 
5  IC1=I 

Xi=X2 

10  CONTINUE 

IF ( IC1) 1 5f 15»  20 
15  RETURN 

20  I  =  I R1 

DO  30  J=1 » NR 
X1T ( J)=E< I> J) 

30  CONTINUE 

I  =  IC1 

DO  40  J=1 » NR 
X2T (J)=E(I>J) 

40  CONTINUE 

I  =  IR1 

DO  50  J=1 » NR 
E< I >  J  > =X2T ( J ) 

50  CONTINUE 

I  =  IC1 

DO  60  J=1 t NR 
E< I > J)=X1T ( J) 

60  CONTINUE 

C 

C  INTERCHANGE  COLUMNS 

C 
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I 


DO  70  1=1 »NR 
XlTd>=Ed,  J> 
CONTINUE 
J=IC1 

DO  80  1=1 » NR 
X2Td>=£d,  J> 
CONTINUE 
J=IR1 

DO  90  1=1 »NR 
Ed,  J)=X2T<I) 
CONTINUE 
J=IC1 

DO  110  1=1, NR 
E(I, J)-XIT(I) 
CONTINUE 
RlT=FdRl> 

R'2T =F dCl ) 
FdRl  > =R2T 
FdCl )  =RlT 
I=IKT(IR1) 
J=IKT  < IC1 ) 

IKT  <  IR1 )  =  J 
IKT dCl  )  =1 
RETURN 
END 
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noon 


PRONY  SUBROUTINE  FOR  CALCULATING  MEAN  SQUARE  DEVIATION 
BETWEEN  PRONY  RECONSTRUCTED  FILE  AND  ACTUAL  DATA  FILE. 

SUBROUTINE  MSO£V( I SP  f  DATA  *  DEV  >  NPTS » ROOTZ  >  ACOEF >  NR ) 
DIMENSION  DAT  A  < 1024 ) 

COMPLEX  ROOTZ ( 15) » ACOEF (15) tB 
C 
C 

DEV=0 . 0 

DO  20  IX=ISP»NPTS+ISP-1 

IT=IX-ISP 

B-CMPLX (0.0»0.0) 

DO  10  IK=1 t NR 

B= ACOEF ( IK ) # ( ROOTZ (IK) **I T  >+B 
10  CONTINUE 

DEV= < REAL ( B > -DATA (IX)) **2+DE V 
20  CONTINUE 

DEV=DEV/NPTS 

RETURN 

END 
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on  nnon  nono 


PRONY  SUBROUTINE  FOR  READING  INPUT  DATA  FILES  IN 
SINGLE  PRECISION  USING  I AG  HEADER  FORMAT 


SUBROUTINE  RFILES<  DATA , J, OX ) 

J=NUMBER  OF  DATA  POINTS  IN  FILE 
DX=TIME  INCREMENT 

DIMENSION  DATA ( 1024) » ICH( 1 ) , IHD (74) 

BYTE  NAME (34) , I TXT (148) » CHAR (2) 

EQUIVALENCE ( CHAR ( 1 ) , ICH < 1 ) ) , ( IHD ( 1 ) *  I TXT ( 1 ) ) 


WRITE ( 6 , 100 ) 

READ < 5» 110) NA» NAME 
NAMECNA+l )=0 
WRITE(2» 120) 

WRITE ( 2  » 130 ) ( NAME (IX) , IX=1 »  NA ) 

OPEN < UNIT=4 , NAME=NAME ,  TYPE* ' OLD ' > FORM* ' UNFORMATTED ' ,  READONLY 
IKC=1 

READ ( 4 , END* 1 SO , ERR* 1 90 ) I , J , K » SX , DX , I CH 

IKC*IKC+1 

WRITE ( 2 » 140 ) 

WRITE < 2  , 150 ) I » J, K ,  SX ,  DX » CHAR( 1 ) >  CHAR ( 2 ) 

NUM=CHAR ( 2 ) 

IF(NUM-1 ) 15 , 15,5 
5  DO  10  IX=1 » NUM+ 1 

READ <4, END* ISO? ERR* 190) IHD (IX) 

10  CONTINUE 

WRITE <2* 160) (ITXT(IX) f IX=1,2*NUM> 

15  IKC*IKC+1 

DO  20  IX=1,J 

READ ( 4 , END= 1 80 , ERR* 1 90 ) DA  TA  < I X ) 

20  CONTINUE 

CALL  CLOSE ( 4 ) 

RETURN 

100  FORMAT ( /  >  '  ENTER  FILE  SPECS.') 

110  FORMAT ( Q ,  34A1 ) 

120  FORMAT ( / » 16X , '  INPUT  DATA  FILE  INFORMATION',/) 

130  FORMAT < 16X, '  ' ,<2*NA>A1 ) 

140  FORMAT ( /, 16X, '  FILE  HEADER  INFORMATION',/) 

150  FORMAT </, 13X, 316 , F12. 5? E12. 5 > 4X , 1A1 ,4X, 14,/) 
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to  to  ro 


160 

180 


FORMAT ( 16X ,  '  '<2*NUM>A1 ) 

WRITE  <  6 ,  200 ) IKC 
GO  TO  220 
190  WRITE(6»210)IKC 

0  FORMAT ( '  END  OF  FILE  ON  READ-PROG.  EXIT", 13,/) 

0  FORMAT ( '  ERR  ON  READ-PROG.  EXIT", 13,/) 

0  CALL  EXIT 

END 


54 


n  n n  no  noon  no 


SUBROUTINE  PROD ( COE , NR ,  I E  R  «  R  0  0  T  Z  ,  I  A)  (see  note  at  end  of 

this  appendix) 

DIMENSIONED  DUMMY  VARIABLES 
DIMENSION  E( 16)  t Q( 16)  » COE (16) t POL <16  > 

COMPLEX  ROOTZ ( 15 ) 

NORMALIZATION  OF  GIVEN  POLYNOMIAL 
TEST  OF  DIMENSION 

IR  CONTAINS  INDEX  OF  HIGHEST  COEFFICIENT 
IER  =  0 

IC=NR-I A+l 
IR  =  IC 
EF‘3  =  1  .  E  -  d 
TOL=i . E-3 
LIMIT=10SIC 
K0UNT=0 

1  IF  < IR- 1 ) 79 , 79  ?  2 


DROP  TRAILING  ZERO  COEFFICIENTS 

2  IF  <  COE ( IR ) )  4  *  3  » 4 

3  IR-IR-i 
GOTO  1 

REARRANGEMENT  OF  GIVEN  POLYNOMIAL 
EXTRACTION  OF  ZERO  ROOTS 

4  0-1 ./COE(IR) 

I E  N  D  =  I R  - 1 
1ST  A=  1 
NSAV=IR+1 
JPEG- 1 

G  C  J )  -•  1  . 

Q  (  J  t  I )  -  C  0  E  ( I R  - 1 )  /  C  0  E  ( I R  ) 

CMirO-COE'  J)/COE(IRj 

WHERE  J  IS  THE  INDEX  OF  THE  LOWEST  NONZERO  COEFFICIENT 

00  9  1=1, IR 

J=NSAV-T 

I F  (COE (1)17,5,7 

5  €0T0<  6,9)  ,  J8EC5 

6  NSAV=NSAV+i 
Q  C  1 5  T  A )  =  0 . 

E(lStA)=0, 

1ST  A-=I  STA  +  1 


55 


n  o 


#  ;o 

^•y 

<?/ 


GOTO  9  ^ 

7  J  6  E  G  =  2  a?  / 

3  Q(  J)=C0E(I'*0  <J  / 

COE I )  =Q  (  J  )  ^  / 

9  CONTINUE  aJT  / 

INITIALIZATION  £  J* 

ESA 0  =  0  . 

Q  < ISTA ) =0  .  £  o' 

10  NSAV=IP  ft  .% 

tJ  ;o 
«>  P 

COMPUTATION  OF  DERIVATIVE  g  .j 

EXPTsIR-ISTA  #  £ 

E  C ISTA ) =EXP T  g  ± 

DO  11  I»  ISTA  fIEND 

EXP  f  =  E:<PT-l  .0  p'  Jr 

PQL<  1  +  1  >=EF'S*ABS(G< 1+1) ) +EPS  *> 

11  E(I+1)=Q(I+1)#EXPT 

TEST  OF  REMAINING  DIMENSION 
I F  ( JSTA-IEND) 12 r  20  r  60 

12  JEND= IEND-1 

COMPUTATION  OF  3-FRACTION 
'DO  1?  I=ISTAfJEND 
IF(I-ISTA>13f16f13 

13  :F(ABS<E(I) )“POL< I  +  l > > 1 4 , 14 f to 

THE  GIVEN  POLYNOMIAL  HAS  MULTIPLE  ROOTS  j  THE  CQEFFICIEN 

THE  COMMON  FACTOR  ARE  STORED  FROM  Q(MSAV)  UP  TO  (KIR) 

14  N3AU-I 

DO  15  K  =  I ?  JEND 

I F ( A B S ( E ( K ) )-PGL(R+l ) ) 1 5  r 1 5  r 8 0 

15  CONTINUE 
GOTO  21 

EUCLIDEAN  ALGORITHM 

16  DO  19  K--I  *  I  END 
ECKrl >*£<K+1 )/E( 1 > 

0  <  f  1.  >  -  E  (  K  +  .1.  >  -  Q  <  K  f  1  ) 

I F  (  K  -  I )  1 3  >  1 7  f  L  8 

TE3T  FOR  SMALL  DIVISOR 

17  IF  C.ABS1QC  1  +  1)1 -POL  ir+L)  >80.80.  19 
IS  GUC+l>aOCK*l>/OU+l) 

POL t K+ 1 ) -POL ( K+l ) /ASS( Q u  + 1 )  ; 
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rj  a  o  o  o 


E ( K ) =G  <  K  +  l > -E ( K ) 

CONTINUE 

Q(IR)=-Q(IR> 

THE  DISPLACEMENT  EXPT  IS  SET  TO  0  AUTOMATICALLY. 

E< ISTA)=0. »Q( ISTA+1 ) , . ♦ . ,E(NSAV-1 >  >G(NSAV> rE(NSAV>=0 

FORM  A  DIAGONAL  OF  THE  QD- ARRAY. 

INITIALIZATION  OF  BOUNDARY  VALUES 
E(ISTA)=0. 

N  R  A  N  =  N  S  A  V  - 1 
E  (  NRANt  i  )  =0  . 


TEST  FOR  LINEAR  OR  CONSTANT  FACTOR 

NRAN-ISTA  IS  DEGREE-1 
I F ( NR AN- 1 STA  >  24  » 23 » 3 1 

LINEAR  FACTOR 
GK  ISTA+1 )=Q< ISTA+1 >+EXPT 
E<  ISTA+1 )=0. 

TEST  FOR  UNFACTORED  COMMON  DIVISOR 
E  < ISTA ) =ESAV 
IF< IR-NSAW)60»A0*25 

INITIALIZE  QD- ALGORITHM  FOR  COMMON  DIVISOR 
I  ST A  =  NS A V 
E  3  A  V  =  E ( ISTA) 

GOTO  10 

COMPUTATION  OF  ROOT  PAIR 
P-P+EXPT 

TEST  FOR  REALITY 
IF ( 0 ) 27 i 28 i 23 

COMPLEX  ROOT  PAIR 
0  NR  AN  >  =P 
Q<NRAM+1  >=P 
E  (  N  R  A  N  >  -  T 
E  (  M  R  A  N  t  1  )  —  “■ 

GOTO 

REAL  ROOT  PAIR 
(3(.Nl?AN)-P-T 
Q(NRAN+1 ) =P+T 
ECNRAN) =0. 


-  Cj 


REDUCTION  OF  DEGREE  BY  2  (DEFLATION) 

2?  N R A N  =  N R A N - 2 
GOTO  22 

COMPUTATION  OF  REAL  ROOT 

3C  Q(NRAN+i )=EXPT+P 

REDUCTION  OF  DEGREE  BY  1  (DEFLATION) 
NRAN=NR AN- 1 
GOTO  22 

START  QD -ITERATION 

31  JBEG=ISTA+1 
J  E  N  D  =  N  R  A  N  - 1 
TEFS=EPS 
T  D  E  L  T  =  1  .  E  -  2 

32  KOUNT=KOUNT  +  1 
P=Q(NRAN+1 ) 

R= ABS ( E ( NRAN ) ) 

TEST  FOR  CONVERGENCE 
IF ( R-TEPS ) 30  f  30  »  33 

33  S-ABS ( E  (  JEND ) ) 

IS  THERE  A  REAL  ROOT  NEXT 
IF ( S-R ) 30  ?  33  »  34 

IS  DISPLACEMENT  SMALL  ENOUGH 

34  IF ( R-TBELT ) 36  r  35  r 35 

35  P*0. 

36  Q-F 

DO  37  J~ JBEG  ?  NRAN 
G(J)-G(J)+E(J>-E<J-1  )~0 

TEST  FOR  SMALL  DIVISOR 
7  F  I A  B  3  (  G  •:  J )  )  -  P  0  L  •:  J  )  >  e  .l .  8  X  f  3 

37  E  (  J  )  ~Q  (  J+l )  *E  (  J  )  /Q  (  J  ) 

0  (  N  R  A  N  +  1  )  =-  -  E  :  N  R  A  N  ’■  -  Q  '  N  R  A  N  *  J.  )  •  -  0 
G  0  T  0  5  4 


C  A C  U  L  A  T  E  j  1 3  P  L  A  C  -I  <1  ENT  F  ( !  R  0  0  !J  B  L  E  R  0  U  7  S 

quadratic  e au at  ion  fir  double  roots 

X  *  *2 -  c  q  (  NR  AN )  +0  ( HR  AN+1  )4€  (  NR  AN)  )  *X+ Q  <  NtfAN )  *G  ( NR  A  N+  1 ; 

38  P=0.5*(Q(W*?^)+E(NRAN)-HJ(NJ?AN  +  3)  ) 

0*P*P-Q(N*AN)#Q(NRAN41) 


o  o  o  o  n  n  n  o 


43  IF ( S-TDELT >44,35,  35 


INITIALIZATION 

44  0=Q(JBEG)+E( JBEG)-P 

TEST  FOR  SMALL  DIVISOR 
IF (ABS(O) -POL  < JBE3 ) ) 3 1 ,  8 1  *  45 

45  T=CT/G>«2 

U  =  EC  JBEG>*CK  JPEG+1  >/<G*(l  ,+T)  ) 

v=o+u 

KOUNT “KOUNT  +  2 

THREEFOLD  LOOP  FOR  COMPLEX  DISPLACEMENT 
DO  53  J= JBEG  ,  NRAN 
0=0 (J+l >+E( J+l ) -U-P 


TEST  FOR  SMALL  DIVISOR 


IF  (AS 

S  (  V  )  -r'OL 

( J 

) 

>46 

,  46  , 

49 

46 

IF(  J- 

H  R  A  N  )  8 1  ? 

4  / 

7 

31 

4? 

E  X  P  T  = 

E  A  P  1  T  P 

I F  (  A  B 

S  (  E  (  J  E  M  D 

) ) 

- 

T'H. 

)  48  ? 

43,31 

i  o 
»  w 

P  -  0  *  5 

T  (  V  +  Q-E  \ 

JE 

a\ 

it 

S  ■  ) 

0  =  P  #  p 

-  (  V  -•  U  )  :!<  ( 

0- 

U 

k  r  - 

Q*U* 

<  1  •  ‘V 

T  =  SC!R 

~ <. A B  5(0 ) 

\ 

SOTO 

L  o 

C  TEST  FOR  SMALL  DIVISOR  J 

49  Jr  (ABS(.0)-PQLCJ+1)  )46,46,:>0 

50  U  =  V*0/l/ 

T=T*(V/0>**2 

<UJ)=V+W-EU~1>  j 

u  -o .  i 


*  ;y  .  . 

Oojri  **/****&&&$ 


r>  n.  n  n  n  nn  n  on  on 


IFCJ-NRAN)51>32?52 

51  U  =  Q( J+2)*E( J  +  l >/<G*< 1  . +T)  ) 

52  V=0+U~W 

TEST  FOR  SMALL  DIVISOR 
IF ( ABS ( CH  J ) ) -POL <  J  > )81»81»53 

53  E( J)=U*V*< 1 ,+T)/Q( J) 

GKNRAN+1  )=V~E(NRAN> 

54  EXPT  =  EXF‘TtP 
TEP3=TEPS*1 .  1 
TDELT=TDELT*1.1 

IF ( KOUNT-LI MIT ) 32 • 55  > 55 

NO  CONVERGENCE  WITH  FEASIBLE  TOLERANCE 

ERROR  RETURN  IN  CASE  OF  UNSATISFACTORY  CONVERGENCE 

55  I E  R  = 1 

REARRANGE  CALCULATED  ROOTS 

56  IENB'-NSAV-NRAN  - 1 
E  ( I S  T  A )  =  E  S  A  V 

I F  < I END) 59 >59 t 57 

57  DO  58  I»lfIEND 
J  ~  I S  T  A + 1 
K=NRAN+1+I 

E( J)=E'!<) 

58  Q ( J ) - Q ( K ) 

59  I R  =  I S T A  + 1 E M D 

NORMAL  RETURN. 

60  I R - I R - 1 

IF  < IR ) 78  f  78  r  61 

REARRANGE  CALCULATED  ROOTS 

61  DO  62  I  - 1 » I R 
GK  I  >  =CK  1+1) 

62  E  I  )  -  E  (  I  +  1  ) 

CALCULATE  COEFFICIENT  VECTOR  '.-ROM  ROOTS 

P0LCIR+1)=1. 

iend=ir-i 

J8E&-1 

DO  69  J*.l  ,  IR 
ISTA=IR+1-J 

0=c , 

PaQl ISTA) 

T*E ( ISTA ) 

IF(T)6S,63 ,65 


trtib 


60 


non  nn  on  o  n 


MULTIPLY  WITH  LINEAR  FACTOR 

63  DO  64  I  =  I S T  A  *  I R 
POL ( I ) =Q-P$POL  ( 1  + 1 ) 

64  0=P0L ( I +1 ) 

GOTO  69 

65  G0T0<66> 67) , JBEG 

66  JBEG=2 

POL ( 1ST  A ) =0  . 

GOTO  69 

MULTIPLY  WITH  QUADRATIC  FACTOR 

67  JBEG= 1 
U=P*P+T*T 
P  =  P  +  P 

DO  68  I  =  I  ST  A 1 1 END 

P0L( I )=0-P*PGL< 1  +  1 )+IJ¥P0L (1+2) 

68  0  =  P0L  ( I  +  .1 ) 

POL ( IR ) “0-P 

69  CONTINUE 

IF ( IER >73,70,78 

COMPARISON  OF  COEFFICIENT  VECTORS i  IE.  TEST  OF  ACCURACY 

70  P=0, 

DO  75  1=1 7 IR 
IF(C0E(I) >72f7l,72 

71  0  =  ABS ( POL ( I )  ) 

GOTO  73 

72  0=ABS< ( POL  < I ) -COE ( I ) ) /COE  ( I  >  ) 

73  IF (P-0) 74. 75,75 

74  P  =  0 

75  CONTINUE 
I F  <  P -- TO L  77 , 76 , 7 6 

76  I ER=~ 1 

77  O’  I R  + 1 ) * P 
£ ( I P+1 )=0 . 

70  DC)  100  I  =  I A  + 1  fN R 

ROOTE ( I )=CMPLX < n< I-IA) •£ ( I-IA ) ) 

100  CONTINUE 
RETURN 

ERROR  RETURNS 

ERROR  RETURN  FOR  POLYNOMIALS  OP  DEGREE  LESS  THAN  1 

79  TER- 2. 

3R=0 
RETURN 


I 


o  n 


ERROR  RETURN  IF  THERE  EXISTS  NO  3-FRACTION 

80  IER=4 
IR= ISTA 
GOTO  60 

ERROR  RETURN  IN  CASE  OF  INSTABLE  QD-ALGORITHM 

81  IER=3 
GOTO  56 
END 


Re  info  on  page  55: 

Taken  from  page  183 
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